%Sensitities
function H_a=Sensitivity(variable,type,f1,g1,f2,g2,x1,y1,z1,p,p1,RM,x,y,z,cn1,cn2,H_Phi,nnod)
sa=variable(:,1);sb=variable(:,2);st=variable(:,3);ca=variable(:,4);cb=variable(:,5);ct=variable(:,6);
x0=variable(:,7);y0=variable(:,8);z0=variable(:,9);L=variable(:,10);a1=variable(:,11);
xp=x1.^p;yp=y1.^p1;zp=z1.^p1;f1p=f1.^(-p1-1);f2p=f2.^(-p1-1);
if type==1      %Circular tube with uniform thickness
    t1=variable(:,12);
    yzp=yp+zp;
    %The sensitivity of Phi with respect to a1
    Phi1_a1=p1*cn1.*f1p.*yzp;
    Phi2_a1=-p1*(1-t1).*cn2.*f2p.*yzp;
    Phi_a1=Phi1_a1+Phi2_a1;
    %The sensitivity of Phi with respect to t1
    Phi_t1=p1*a1.*cn2.*f2p.*yzp;
    %The sensitivities of f1,g1,f2 and g2 with respect to x1
    f1_x1=sparse(nnod,1); g1_x1=sparse(nnod,1);
    f2_x1=sparse(nnod,1); g2_x1=sparse(nnod,1);
elseif type==2   %Circular tube with linearly varying thickness
    t1=variable(:,12);t2=variable(:,13);
    yzp=yp+zp;
    %The sensitivity of Phi with respect to a1
    Phi1_a1=p1*cn1.*f1p.*yzp;
    Phi2_a1=-p1*cn2.*f2p.*yzp./a1;
    Phi_a1=Phi1_a1+Phi2_a1;
    %The sensitivity of Phi with respect to t1
    Phi_t1=p1*cn2.*a1.*f2p*(-0.5+0.5*x1./L).*yzp;
    %The sensitivity of Phi with respect to t2
    Phi_t2=p1*cn2.*a1.*f2p*(-0.5-0.5*x1./L).*yzp;
    %The sensitivities of f1,g1,f2 and g2 with respect to x1
    f1_x1=sparse(nnod,1); g1_x1=sparse(nnod,1);
    f2_x1=0.5*(t1-t2)./L.*a1; g2_x1=f2_x1;
elseif type==3    %Tapered circle tube with uniform thickness
    a2=variable(:,12);t1=variable(:,13);
    yzp=yp+zp;
    %The sensitivities of Phi with respect to a1 and a2
    Phi1_a1=p1*cn1.*f1p.*yzp.*((1+a2)-(1-a2).*x1./L)*0.5;
    Phi1_a2=p1*cn1.*a1*f1p.*yzp.*(1+x1./L)*0.5;
    Phi2_a1=-p1*cn2.*f2p.*yzp.*((1+a2)-(1-a2).*x1./L-2*t1.*a2)*0.5;
    Phi2_a2=-p1*cn2.*a1*f2p.*yzp.*(1+x1./L-2*t1)*0.5;
    Phi_a1=Phi1_a1+Phi2_a1;
    Phi_a2=Phi1_a2+Phi2_a2;
    %The sensitivity of Phi with respect to t1
    Phi_t1=p1*cn2.*a1*a2*f2p*yzp;
    %The sensitivities of f1,g1,f2 and g2 with respect to x1
    f1_x1=0.5*(a1.*a2-a1)./L;g1_x1=f1_x1;
    f2_x1=f1_x1;g2_x1=f1_x1;
    
elseif type==4 %Rectangular tube with uniform thickness
    b1=variable(:,12);t1=variable(:,13);
    g1p=g1.^(-p-1);   g2p=g2.^(-p-1);
    %The sensitivities of Phi with respect to a1 and b1
    Phi1_a1=p*cn1.*(f1p.*yp+b1.*g1p.*zp);
    Phi1_b1=p*cn1.*a1.*g1p.*zp;
    Phi2_a1=-p*cn2.*((1-t1.*b1).*f2p.*yp-(t1.*b1-b1).*g2p.*zp);
    Phi2_b1=-p*cn2.*((a1-t1.*a1).*g2p.*zp-t1.*a1.*f2p.*yp);
    Phi_a1=Phi1_a1+Phi2_a1;
    Phi_b1=Phi1_b1+Phi2_b1;
    %The sensitivity of Phi with respect to t1
    Phi_t1=p*cn2.*a1.*b1.*(f2p.*yp+g2p.*zp);
    %The sensitivities of f1,g1,f2 and g2 with respect to x1
    f1_x1=sparse(nnod,1); g1_x1=sparse(nnod,1);
    f2_x1=sparse(nnod,1); g2_x1=sparse(nnod,1);
    
elseif type==5  %Rectangular tube with linearly varying thickness
    b1=variable(:,12);t1=variable(:,13);t2=variable(:,14);
    g1p=g1.^(-p-1);g2p=g2.^(-p-1);
    cot=a1.*b1; co=0.5*(t1+t2)+0.5*(t2-t1)./L;
    %The sensitivities of Phi with respect to a1 and b1
    Phi1_a1=p*cn1.*(f1p.*yp+b1.*g1p.*zp);
    Phi1_b1=p*cn1.*(g1p.*zp+a1.*f1p.*yp);
    Phi2_a1=p*cn2.*((b1.*co-1).*b1.*f2p.*yp+co.*zp.*g2p);
    Phi2_b1=p*cn2.*((a1.*co-1).*a1.*g2p.*zp+co.*yp.*f2p);
    Phi_a1=Phi1_a1+Phi2_a1;
    Phi_b1=Phi1_b1+Phi2_b1;
    %The sensitivities of Phi with respect to t1 and t2
    Phi_t1=p*cn2.*cot.*(1/2-x1/2./L)*(f2p.*yp+g2p.*zp);
    Phi_t2=p*cn2.*cot.*(1/2+x1/2./L)*(f2p.*yp+g2p.*zp);
    %The sensitivities of f1,g1,f2 and g2 with respect to x1
    f1_x1=sparse(nnod,1); g1_x1=sparse(nnod,1);
    f2_x1=-(t2-t1)/2./L.*cot; g2_x1=f2_x1;
    
elseif type==6 %%Tapered rectangular tube with uniform thickness
    a2=variable(:,12);b1=variable(:,13);b2=variable(:,14);t1=variable(:,15);
    co1=0.5-0.5*x1./L;                 co2=0.5+0.5*x1./L;
    co11=0.5*(1+a2)-0.5*(1-a2).*x1./L; co22=0.5*(b1+b2*a2)-0.5*(b1-b2.*a2).*x1./L;
    g1p=g1.^(-p-1);g2p=g2.^(-p-1);
    %The sensitivities of Phi with respect to a1, a2, b1 and b2
    Phi1_a1=p*cn1.*(f1p.*yp.*co11+g1p.*zp.*co22);
    Phi1_a2=p*cn1.*(a1.*f1p.*yp.*co2+a1.*b2.*g1p.*zp.*co2);
    Phi1_b1=p*cn1.*a1.*g1p.*zp.*co1;
    Phi1_b2=p*cn1.*a1.*a2.*g1p.*zp.*co2;
    Phi2_a1=-p*cn2.*(f2p.*yp.*(co11-a2.*b2.*t1)-g2p.*zp.*(co22-a2.*b2.*t1));
    Phi2_a2=-p*cn2.*(a1.*f2p.*yp.*(co2-b2.*t1)-a1.*b2.*g2p.*zp.*(co2-t1));
    Phi2_b1=-p*cn2.*a1.*g2p.*zp.*co1;
    Phi2_b2=-p*cn2.*a1.*a2.*(g2p.*zp.*(co2-t1)-t1.*f2p.*yp);
    Phi_a1=Phi1_a1+Phi2_a1; Phi_a2=Phi1_a2+Phi2_a2;
    Phi_b1=Phi1_b1+Phi2_b1; Phi_b2=Phi1_b2+Phi2_b2;
    %The sensitivity of Phi with respect to t1
    Phi_t1=p*a1.*b1.*(f2p.*yp+g2p.*zp);
    %The sensitivities of f1,g1,f2 and g2 with respect to x1
    f1_x1=0.5*(a1.*a2-a1)./L;dg1_x1=0.5*(a1.*a2.*b2-a1.*b1)./L;
    f2_x1=f1_x1;g2_x1=dg1_x1;
end

%The sensitivity of Phi with respect to L, alfa, beta, theta, x0, y0 and z0
xx0=x-x0;yy0=y-y0;zz0=z-z0;
bjt1=p*(x1./L).^(p-1)./L;
bjt2=p1*(y1./f1).^(p1-1)./f1./f1;
bjt3=p1*(z1./g1).^(p1-1)./g1./g1;
bjt4=p1*(y1./f2).^(p1-1)./f2./f2;
bjt5=p1*(z1./g2).^(p1-1)./g2./g2;
dT_alfa{1}=0;dT_alfa{2}=0;dT_alfa{3}=0;
dT_alfa{4}=ca.*sb.*ct-sa.*st; dT_alfa{5}=-ca.*sb.*st-sa.*ct; dT_alfa{6}=-ca.*cb;
dT_alfa{7}=sa.*sb.*ct+sa.*st; dT_alfa{8}=-sa.*sb.*st+sa.*ct; dT_alfa{9}=-sa.*cb;
dT_beta{1}= -sb.*ct;dT_beta{2}=sb.*st ;dT_beta{3}=cb ;
dT_beta{4}= sa.*cb.*ct;dT_beta{5}=-sa.*cb.*st ;dT_beta{6}= sa.*sb;
dT_beta{7}= -ca.*cb.*ct;dT_beta{8}= ca.*cb.*st;dT_beta{9}= -ca.*sb ;
dT_theta{1}=-cb.*st;dT_theta{2}=-cb.*ct;dT_theta{3}=0;
dT_theta{4}=-sa.*sb.*st+ca.*ct ;dT_theta{5}=-sa.*sb.*ct-ca.*st ;dT_theta{6}=0;
dT_theta{7}=ca.*sb.*st+sa.*ct ;dT_theta{8}=ca.*sb.*ct-sa.*st ;dT_theta{9}=0;
dx1_alfa=xx0.*dT_alfa{1}+yy0.*dT_alfa{2}+zz0.*dT_alfa{3};
dy1_alfa=xx0.*dT_alfa{4}+yy0.*dT_alfa{5}+zz0.*dT_alfa{6};
dz1_alfa=xx0.*dT_alfa{7}+yy0.*dT_alfa{8}+zz0.*dT_alfa{9};
dx1_beta=xx0.*dT_beta{1}+yy0.*dT_beta{2}+zz0.*dT_beta{3};
dy1_beta=xx0.*dT_beta{4}+yy0.*dT_beta{5}+zz0.*dT_beta{6};
dz1_beta=xx0.*dT_beta{7}+yy0.*dT_beta{8}+zz0.*dT_beta{9};
dx1_theta=xx0.*dT_theta{1}+yy0.*dT_theta{2}+zz0.*dT_theta{3};
dy1_theta=xx0.*dT_theta{4}+yy0.*dT_theta{5}+zz0.*dT_theta{6};
dz1_theta=xx0.*dT_theta{7}+yy0.*dT_theta{8}+zz0.*dT_theta{9};
dx1_x0=-RM{1};dy1_x0=-RM{4};dz1_x0=-RM{7};
dx1_y0=-RM{2};dy1_y0=-RM{5};dz1_y0=-RM{8};
dx1_z0=-RM{3};dy1_z0=-RM{6};dz1_z0=-RM{9};
dPhi1_L= p*cn1.*xp.*L.^(-p-1);
dPhi1_alfa=cn1.*(-bjt1.*dx1_alfa-bjt2.*(dy1_alfa.*f1-f1_x1.*dx1_alfa.*y1)-bjt3.*(dz1_alfa.*g1-g1_x1.*dx1_alfa.*z1));
dPhi1_beta=cn1.*(-bjt1.*dx1_beta-bjt2.*(dy1_beta.*f1-f1_x1.*dx1_beta.*y1)-bjt3.*(dz1_beta.*g1-g1_x1.*dx1_beta.*z1));
dPhi1_theta=cn1.*(-bjt1.*dx1_theta-bjt2.*(dy1_theta.*f1-f1_x1.*dx1_theta.*y1)-bjt3.*(dz1_theta.*g1-g1_x1.*dx1_theta.*z1));
dPhi1_x0=cn1.*(-bjt1.*dx1_x0-bjt2.*(dy1_x0.*f1-f1_x1.*dx1_x0.*y1)-bjt3.*(dz1_x0.*g1-g1_x1.*dx1_x0.*z1));
dPhi1_y0=cn1.*(-bjt1.*dx1_y0-bjt2.*(dy1_y0.*f1-f1_x1.*dx1_y0.*y1)-bjt3.*(dz1_y0.*g1-g1_x1.*dx1_y0.*z1));
dPhi1_z0=cn1.*(-bjt1.*dx1_z0-bjt2.*(dy1_z0.*f1-f1_x1.*dx1_z0.*y1)-bjt3.*(dz1_z0.*g1-g1_x1.*dx1_z0.*z1));
dPhi2_L= -p*cn2.*xp.*L.^(-p-1);
dPhi2_alfa=cn2.*(bjt1.*dx1_alfa+bjt4.*(dy1_alfa.*f2-f2_x1.*dx1_alfa.*y1)+bjt5.*(dz1_alfa.*g2-g2_x1.*dx1_alfa.*z1));
dPhi2_beta=cn2.*(bjt1.*dx1_beta+bjt4.*(dy1_beta.*f2-f2_x1.*dx1_beta.*y1)+bjt5.*(dz1_beta.*g2-g2_x1.*dx1_beta.*z1));
dPhi2_theta=cn2.*(bjt1.*dx1_theta+bjt4.*(dy1_theta.*f2-f2_x1.*dx1_theta.*y1)+bjt5.*(dz1_theta.*g2-g2_x1.*dx1_theta.*z1));
dPhi2_x0=cn2.*(bjt1.*dx1_x0+bjt4.*(dy1_x0.*f2-f2_x1.*dx1_x0.*y1)+bjt5.*(dz1_x0.*g2-g2_x1.*dx1_x0.*z1));
dPhi2_y0=cn2.*(bjt1.*dx1_y0+bjt4.*(dy1_y0.*f2-f2_x1.*dx1_y0.*y1)+bjt5.*(dz1_y0.*g2-g2_x1.*dx1_y0.*z1));
dPhi2_z0=cn2.*(bjt1.*dx1_z0+bjt4.*(dy1_z0.*f2-f2_x1.*dx1_z0.*y1)+bjt5.*(dz1_z0.*g2-g2_x1.*dx1_z0.*z1));
dPhi_L=dPhi1_L+dPhi2_L;
dPhi_alfa=dPhi1_alfa+dPhi2_alfa; dPhi_beta=dPhi1_beta+dPhi2_beta; dPhi_theta=dPhi1_theta+dPhi2_theta;
dPhi_x0=dPhi1_x0+dPhi2_x0;       dPhi_y0=dPhi1_y0+dPhi2_y0;       dPhi_z0=dPhi1_z0+dPhi2_z0;
%The sensitivity of H with respect to design variables
H_a{1}= H_Phi.*dPhi_alfa;H_a{2}= H_Phi.*dPhi_beta;H_a{3}= H_Phi.*dPhi_theta;
H_a{4}= H_Phi.*dPhi_x0;  H_a{5}= H_Phi.*dPhi_y0;  H_a{6}= H_Phi.*dPhi_z0;
H_a{7}= H_Phi.*dPhi_L;   H_a{8}= H_Phi.*Phi_a1;
if type==1
    H_a{9}= H_Phi.*Phi_t1;
elseif type==2
    H_a{9}= H_Phi.*Phi_t1;H_a{10}=H_Phi.*Phi_t2;
elseif type==3
    H_a{9}=H_Phi.*Phi_a2; H_a{10}=H_Phi.*Phi_t1;
elseif type==4
    H_a{9}=H_Phi.*Phi_b1; H_a{10}=H_Phi.*Phi_t1;
elseif type==5
    H_a{9}=H_Phi.*Phi_b1; H_a{10}=H_Phi.*Phi_t1; H_a{11}=H_Phi.*Phi_t2;
elseif type==6
    H_a{9}=H_Phi.*Phi_a2; H_a{10}=H_Phi.*Phi_b1; H_a{11}=H_Phi.*Phi_b2;  H_a{12}=H_Phi.*Phi_t1;
end
end


